Frontiers in Systems Biology

Mechanistic mathematical models have been used extensively to provide a deeper understanding of biological mechanisms, including unveiling the regulation of tumour growth and its response to various treatments. However, given the breadth of biological regulatory mechanisms, these models are frequently large and thus prone to potential issues with parameter identifiability. Statistical metrics like the Akaike and Bayesian information criteria can help identify a parsimonious model by balancing goodness of fit against model complexity. Yet simple models may fail to provide sufficient biological insight if they do not adequately capture known physiological processes or mechanisms. A modeller must therefore balance hypothesis generation and biological learning with model tractability. Here, we illustrate this balance using models of ovarian cancer growth and treatment response to cisplatin and immune checkpoint blockade in homologous recombination (HR)-deficient and HR-proficient immunocompetent mouse models. We develop a hierarchy of mathematical models of increasing complexity to describe tumour growth, treatment response, and immune dynamics. Our results highlight the limits of relying purely on statistical metrics for model selection, particularly when the goal is to obtain biological insight and underscore the importance of balancing model complexity to avoid overfitting and parameter unidentifiability.

Analysis of circadian patterns in time-series data requires computational methods that can accommodate several factors, including variable sampling resolution, replicate number, and missing values. Most existing tools simplify rhythmicity to a strict dichotomy based solely on a single p-value threshold. This leads to a level of uncertainty that affects many biological targets. We developed CODAC (Circadian Oscillation Detection Analysis and Comparison), a framework that integrates nonlinear constrained optimization with a multicriteria rhythmicity classification scheme to evaluate rhythmic patterns without relying on a single statistical cutoff. This approach allows CODAC to identify and exclude medium-confidence rhythms rather than force them into a rhythmic/arrhythmic dichotomy. CODAC comprises four modules: (i) CODAC_single estimates rhythmicity within a single group; (ii) CODAC_flex extends this to identify distinct waveform types within one group; (iii) CODAC_compare performs pairwise comparisons across two or more groups to detect rhythmic or arrhythmic changes; and (iv) CODAC_multi handles more complex designs involving multiple-group comparisons. Using in silico simulations and public transcriptomic datasets, we show that CODAC performs comparably to established methods while providing additional flexibility for rhythm classification and comparison. Taken together, CODAC provides a flexible and open-source package for circadian timeseries analysis with automated visualization tools.

Protein-protein interactions (PPIs) play a crucial role in enabling proteins to carry out their functions within various biological processes (Hui et al., 2003). Since the introduction of the yeast two-hybrid (Y2H) method for PPI detection in 1989 (Fields and Song, 1989), the identification of PPIs has become a significant focus in modern biological research. PPI goes beyond examining individual proteins, allowing researchers to establish a comprehensive network that regulates biological processes. Rice, as a key model organism in plant biological studies, has been at the forefront of PPI research. In 2008, prominent rice scientists in China called for concerted efforts to define a comprehensive protein-protein interaction network experimentally, which aimed to facilitate the prediction of the functional mechanisms operating throughout a plants lifecycle (Zhang et al., 2008). With efforts for 2 decades, the experimentally identified rice PPIs have reached over ten thousand. Several public databases have been established to systematically collate and store PPIs, including STRING (Szklarczyk et al., 2019), BioGRID (Oughtred et al., 2020), IntAct (del Toro et al., 2022), PRIN (Gu et al., 2011), RicePPINet (Liu et al., 2017) and RiceNet v2 (Lee et al., 2015). However, most PPI datasets in rice stem from computational predictions, while experiment-based rice PPI datasets are fragmented due to the lack of systematic profiling at the rice PPIome level, which largely hinders information sharing in the rice research community. To bridge this gap, we constructed the Port of Protein-Protein Interactomes (POPPIN; https://riceome.hzau.edu.cn/poppin/), an integrated database dedicated to sharing experimentally verified PPIs and functional clues in rice. Empowered by high-throughput PPIome profiling technologies and text mining assisted by a large language model (Huang et al., 2025; Liu et al., 2025), POPPIN currently has deposited over 150,451 pieces of rice PPI-related information. Additionally, POPPIN provides detailed protein information, including GO annotations, subcellular localizations, domains, trait ontology (TO) information, and hyperlinks to external biological databases. Through offering a user-friendly web interface for search and dynamic network visualization, POPPIN serves as the first large-scale, experiment-based database for searchable PPIs in rice, and has the potential to be extended to other species under this structural framework.

Rule-based languages have proven effective for modeling systems of interacting structured entities as typically encountered in chemistry and molecular biology. We present PyKappa, a rule-based modeling package written in Python whose interpreted nature enables interactive simulation and analysis, including by agentic AI. The package seeks to broaden the base of developers by utilizing a widely known programming language and serves as an easy-to-deploy teaching tool. Using PyKappa, we conduct a case study of phase separation.

Mechanistic models of photosynthesis have expanded substantially over the past decades, covering processes from light reactions to carbon fixation. However, these models remain fragmented across the literature, inconsistently implemented, and difficult to reproduce or reuse, limiting their adoption beyond the research group that developed them. Here, we present GreenSloth, a freely accessible web-based database of 22 published mechanistic photosynthesis models, reimplemented as standardized, executable Python objects within MxlPy, an open-source framework for mechanistic biological modeling. Although the database is primarily designed for dynamic mechanistic models formulated as ordinary differential equations, the current implementation also includes the fields most widely cited steady-state mechanistic model and its variants. GreenSloth provides a structured environment for model discovery, comparison, and reuse, addressing reproducibility challenges in the field and enabling integration into emerging hybrid modeling approaches. It is also interactive: each model runs directly in the browser, with no installation, environment setup, or programming required. The resource is openly accessible and designed for long-term community maintenance, hoping to position itself as foundational infrastructure for the photosynthesis modeling community. Database URLhttps://greensloth.rwth-aachen.de/

Muscle regeneration relies on the coordinated activation of muscle stem cells, whose fate decisions are regulated by intracellular gene expression dynamics and intercellular coupling via the Notch-Dll1 signaling pathway. Central components of this regulatory network include the transcriptional repressor Hes1, its target gene Dll1, and the myogenic regulator MyoD. Experimental and theoretical studies have shown that proliferating muscle stem cells exhibit oscillatory dynamics of these molecules, whereas sustained expression is associated with differentiation. Here, we investigate the dynamics of a previously established delay differential equation model of two coupled muscle stem cells. Using linear stability analysis, we systematically characterize how model parameters affect the transition between stable and unstable steady states. In addition, numerical bifurcation analysis is employed to study the influence of intercellular coupling strength and delay on the system dynamics. Our analysis shows that continuous variation of the coupling delay induces repetitive changes in the stability of the steady state. However, this sensitivity towards the coupling delay is confined to a narrow region of parameter space and therefore requires a fine tuning of all other parameters. Beside the identification of parameter sets for in-phase and out-of-phase oscillations, we demonstrate the possibility of coexisting stable in-phase and out-of-phase oscillations, a dynamical feature that has not been reported previously. While oscillation periods are largely determined by intracellular regulatory mechanisms, oscillation amplitudes can be strongly modulated by intercellular coupling. These results provide new insight into how the intracellular network and intercellular communication interact to generate different collective dynamics.

Oscillatory and saturating dynamics in regulatory gene networks govern processes such as cell cycle timing, circadian rhythms, and inflammatory signaling, and sparse regression methods such as SINDy offer an interpretable method for recovering their governing equations from time-series data. However, these networks are typically governed by Hill-type nonlinearities that become nearly collinear with polynomial candidates over most bounds of an oscillating curve, and placing them in a shared library produces ill-conditioned regressions and unreliable coefficient estimates. Standard approaches such as ridge regularization, library normalization, and orthogonalization can marginally improve conditioning but cannot remove the structural multicollinearity of these systems, in turn limiting our ability to produce accurate parameter estimates. To address this issue, we propose the Regulation-Informed Sparse Identification of Nonlinear Dynamics (RI-SINDy) method, which removes known regulatory terms from the regression and recovers their coefficients directly from an asymptotic drift-balance condition. Because the time averaged derivative of a bounded trajectory approaches zero, the coefficients of regulatory nonlinearities can be computed independently by constraining them to offset the mean contribution of the remaining terms. We demonstrate that this method effectively reduces library collinearity and improves conditioning numbers by 3-45x, substantially improving coefficient recovery and regression accuracy. We test the RI-SINDy method on experimental Hes1 data and on synthetic Goodwin and NF-{kappa}B models. On Hes1, it achieves the lowest trajectory error among the compared methods (Traditional SINDy, SR3, and Nullcline Reconstruction). On the synthetic systems, it recovers the correct sparse structure with lower trajectory error than SR3 and Nullcline Reconstruction, and in a Goodwin noise sweep it recovers all active terms at noise levels up to 15%, where Traditional SINDy and E-SINDy fail to. In conclusion, the RI-SINDy method provides a robust framework to overcome the issue of multicollinearity for system identification in biological systems.

Selecting stochastic gene-expression models from single-cell counts requires accurate parameter inference and efficient model selection. Likelihood methods in count space can be costly when full stationary count distributions are unavailable, whereas approximate methods may lose accuracy. Probability generating functions (PGFs) offer a compact analytical alternative, but existing PGF workflows are generally not likelihood based and therefore rely on computationally intensive cross-validation. We develop a likelihood-based PGF framework for both tasks. Correlated empirical PGF values are used to construct a Gaussian quasi-likelihood for parameter inference and PGF-based Bayesian information criterion (BIC) for model selection. We show that the empirical PGF is exactly unbiased and that the parameter estimator is consistent, converges at the inverse-square-root sample-size rate, and is first-order asymptotically unbiased. For large samples and a uniquely preferred model, PGF-BIC selects the same model as leave-one-out cross-validation in PGF space.

Biochemical Systems Theory (BST) represents nonlinear biochemical rate laws by local power-law approximations in logarithmic concentration coordinates. First-order coefficients are elasticities, whereas higher-order derivatives describe local log-synergism and its variation. Ordinary higher derivatives, however, are not tensorial under nonlinear reparameterizations and can mix biochemical response structure with coordinate artifacts. We formulate a covariant hierarchy, cBST1-cBST3, on a positive operating-point space equipped with a declared reference connection. cBST1 recovers classical elasticities in a flat logarithmic chart, cBST2 is the covariant Hessian of the log-response, and cBST3 is the symmetrized covariant derivative of cBST2. The framework is quantitatively evaluated using three representative rate laws from the curated yeast glycolysis model BIOMD0000000064: glucose transport, glucose phosphorylation, and phosphofructokinase. For 10,000 finite log-concentration perturbations at each of five radii, cBST2 reduced the cBST1 log-rate root-mean-square error by 96.6-98.6% at the largest tested radius, and cBST3 provided a further 68.4-98.1% reduction. The contracted cBST2 and cBST3 terms strongly predicted the corresponding lower-order truncation errors. Under the nonlinear transformation qi = sinh(ui), covariant contractions agreed across coordinates to within a 95th-percentile relative error of 1.2 x 10-13, whereas ordinary higher derivatives showed order-unity coordinate mismatches. Supplementary analytic tests recovered the expected second-, third-, and fourth-order truncation-error scaling. These results show that cBST1-cBST3 are not only coordinate-consistent descriptors but also practical diagnostics of where local power-law approximations require higher-order correction.

Spatial transcriptomics keeps measurement of gene expression while preserving spatial context, yet traditional analysis methods face challenges in computational efficiency, biological interpretability, and autonomous discovery. This project presents a framework solving these issues through three parts: (1) an ensemble clustering system achieving 66.7% improvement over baseline average and 23.9% over best single method with silhouette score of 0.540 and statistical significance (p = 0.0032, Cohens d = 1.82); (2) a knowledge-based clustering framework that annotates 88.6% of cells across 8 ovarian cell types using 428 marker genes; and (3) a GPT-4o-mini-powered autonomous agent that generated 3 biological hypotheses with validations.

Dose-regimen design requires choosing how much drug to give, when to give it, and how treatment should vary across patients. Mechanistic pharmacokinetic-pharmacodynamic (PK/PD) and quantitative systems pharmacology (QSP) models can predict treatment responses, but optimizing dosing inputs depends on model-specific sensitivity derivations or derivative-free search. Here, we introduce DiffDose, a differentiable programming framework for mechanistic open-loop dose-regimen optimization that uses automatic differentiation (AD) to handle clinically interpretable dose amounts and administration times as differentiable controls. We evaluate our method in three settings: fixed-schedule dose-amplitude optimization in OptiDose PK/PD benchmarks; dose-timing optimization in a chemotherapy-induced neutropenia model with state-dependent delay; and individualized mosunetuzumab dosing in a QSP virtual population. Across these examples, AD produced gradients consistent with references, reduced model-specific derivative work, and shortened benchmark time to solution. DiffDose thereby turns mechanistic PK/PD and QSP models from tools that evaluate prespecified regimens into gradient-based engines for individualized regimen design.

Long-term in vivo transcriptomic time courses are costly, limiting assessment of chronic molecular responses from short studies. We developed a pharmacokinetics-informed transcriptomic ordinary differential equation model (PKT-ODE) that links an oral pharmacokinetic profile and Hill drug-effect function to first-order turnover of co-expression modules. The model was fitted to rat liver responses to fenofibrate at three doses in Open TG-GATEs through day 8. At the held-out day-29 endpoint, PKT-ODE achieved Pearson r = 0.960 and mean squared error (MSE) = 0.148. In this dataset, these values achieved lower prediction error and higher correlation than four statistical baselines and validation-selected linear and multilayer-perceptron transition models. Literature-curated peroxisome proliferator-activated receptor target genes occurred only in modules with positive fitted drug effects. These results provide a proof of concept for pharmacokinetics-informed transcriptomic extrapolation; cross-compound, cross-organ and alternative-regimen performance remain to be tested.

Under current systemic treatment of metastatic cancer, a drug is frequently prescribed at maximum tolerable dose (MTD) until either unacceptable toxicity or progression. Unfortunately, in many patients this treatment strategy leads to the development of treatment resistance. Evolutionary therapy approaches aim to forestall or delay treatment resistance in cancer by exploiting eco-evolutionary interactions. A well-known implementation is the adaptive therapy protocol of Zhang et al., in which tumour burden thresholds are used to guide strategic treatment holidays. Deep reinforcement learning (DRL) has recently been used to optimise these approaches. However, research combining DRL with evolutionary therapy approaches has so far focused on time to progression (TTP) as a performance metric, and has not quantified safety in terms of robustness to delayed treatment restart or included patient preferences regarding quality of life (QoL) in treatment design. In our study, we use a DRL agent informed by a mathematical two-population tumour growth model to design treatment schedules for patients with non-small cell lung cancer (NSCLC). The agent is trained on a virtual patient cohort using parameters previously fitted to data from patients with NSCLC treated with erlotinib. Beyond TTP, we focus on improving robustness to delayed treatment restart and on how individual preferences and values impact QoL experienced during treatment. We compare TTP, robustness and QoL under the DRL policy, the adaptive therapy protocol of Zhang et al., and MTD. We introduce a robustness metric ``margin-to-failure'' (MTF), and compare quality-adjusted-survival (QAS) across different patient preference profiles. Finally, we explore reward shaping to assess how QoL preferences can be incorporated into DRL-based treatment design. To evaluate our results, we consider different decision intervals, defined as the time between dosing adjustments. The DRL policy achieved greater median TTP, MTF, and QAS across all treatment decision intervals compared to the other two protocols. As decision intervals increased, TTP under DRL declined gradually towards that achieved under MTD. In contrast, the Zhang et al. protocol performed inconsistently and could result in premature progression. Additionally, a population-level policy trained on a cohort of virtual patients produced an interpretable treatment rule that extended TTP for most previously unseen patients and indicated that treatment should resume at a lower tumour burden when monitoring is less frequent. These findings show that DRL can balance the benefit of preserving drug-sensitive cells to suppress resistance against the risk of unsafe tumour regrowth. Reward shaping further showed how treatment strategies could be adjusted to reflect different patient preferences. Together, these results provide a biologically informed approach for designing robust and patient-centred evolutionary therapies in fast-growing cancers such as NSCLC.

Metastasis involves context-dependent molecular interactions in which non-coding RNAs, particularly miRNAs and circRNAs, play important regulatory roles. However, existing computational approaches generally do not jointly represent cancer type, metastatic event, and cancer-specific metastatic context. We developed a context-aware multi-task heterogeneous graph neural network (GNN) for predicting ncRNA associations with cancer types and metastatic events. The framework integrates multiple biological repositories into a heterogeneous graph representing ncRNAs, cancers, metastatic event types (METs), and cancer-specific metastatic instances (CSMIs). The model performs six link-prediction tasks using a hierarchical transformer-based encoder and multi-relational TuckER decoder. Across ten independently initialized runs evaluated on the RNA-group-disjoint held-out test set, the model achieved a global AUROC of 0.8801 {+/-} 0.0118 and an F1 score of 0.8260 {+/-} 0.0071. All three ablation variants yielded lower AUROC, with the largest reduction under independent task training. Case studies in pancreatic cancer, colorectal cancer, and hepatocellular carcinoma provided disease-level, event-level, and expression-based support, respectively, for top-ranked candidate associations. The framework enables context-specific prioritization of ncRNA-cancer-metastasis associations for experimental evaluation.

Light history leaves persistent changes in circadian period, but where this history is stored remains unresolved. Suprachiasmatic nucleus (SCN) network models have often approached photoperiodic encoding through phase organization or coupling strength. We computationally tested slow adaptation of subregion-specific intrinsic periods as an alternative memory mechanism. The model asymmetrically couples dorsal (D) and ventral (V) SCN oscillators and adds a systemic oscillator (X) representing putative circadian feedback present in vivo but lost ex vivo. With a single parameter set, period adaptation captured the direction and approximate magnitude of behavioral aftereffects across photoperiod and T-cycle conditions. Adapting coupling strength instead of period failed to reproduce the V-leading-D phase order reported after T22. Removing systemic feedback preserved the photoperiod-dependent period ordering but inverted the T22 and T26 aftereffects, an inversion that matched SCN explant observations. The model also yielded distinct D-V phase organization for each of 18:6 LD, T23, and T25. These results suggest that subregion-specific period plasticity provides a parsimonious substrate for encoding light history, while the dependence on systemic feedback indicates that behavioral period may not be a readout of the SCN alone. Graphical Abstract O_FIG O_LINKSMALLFIG WIDTH=200 HEIGHT=160 SRC="FIGDIR/small/743784v1_ufig1.gif" ALT="Figure 1"> View larger version (33K): org.highwire.dtl.DTLVardef@d48266org.highwire.dtl.DTLVardef@1bd15bdorg.highwire.dtl.DTLVardef@de3d9forg.highwire.dtl.DTLVardef@9fbc39_HPS_FORMAT_FIGEXP M_FIG C_FIG A model with dorsal period adaptation and phenomenological systemic feedback accounts for behavioral and explanted SCN aftereffects. (A) During T22 entrainment, dorsal (D), ventral (V), and systemic (X) oscillators remain phase-locked. After release into constant darkness, systemic coupling maintains a unified in vivo rhythm, whereas removing X feedback in the explant simulation allows the D-V network to express a distinct ex vivo period aftereffect. (B) The SCN model is modeled as an asymmetrically coupled attractive-repulsive oscillator network with stronger photic input to V. Light history is encoded via plasticity of the intrinsic period in D, while X represents putative systemic circadian feedback available in vivo and lacking direct photic input. (C) The model reproduces concordant period changes in behavior and SCN explants across photoperiods, but opposing period changes following T-cycle entrainment.

Scientific computing has become a central component of modern scientific discovery. Yet many computational tools are developed by small, specialized teams under incentives that encourage the release of rapidly prototyped tooling without commensurate attention to engineering concerns, including performance and maintainability. These gaps are particularly visible in the life sciences, where the advent of high-throughput sequencing and molecular profiling has made the production and processing of datasets routine at scales that strain reliability and cost. Recently, LLM-based agents have become increasingly capable, with publicly available systems possessing both significant domain knowledge in many scientific fields and the ability to autonomously operate over complex and specialized codebases in pursuit of well-defined goals. Together, these developments create a practical opportunity for scientific computing. Many of the persistent weaknesses of the scientific computing ecosystem stem from technical debt and a shortage of sustained engineering labor and expertise. Here, we examine coding agents as a potential way to address these weaknesses: we present an exploratory field report of eight early case studies in the application of LLM agents to scientific computing across a range of project scopes, from lightweight maintenance tasks to full performance-oriented rewrites of scientific libraries, with a focus on the life sciences. Each of these case studies is accompanied by reflections from the individual or group responsible for the work, including lessons from the process. Overall, we find that the use of coding agents in scientific computing holds great promise for accelerating scientific research and increasing the reliability of critical systems, but that outstanding concerns remain, including responsibility and ownership for such projects, and we suggest collaboration and stewardship with existing maintainers when feasible.

Cross-dataset generalization of EEG-based classification under weak, proxy-derived labels remains an open problem for altered-states research. We present a reproducible eight-dataset alignment pipeline that maps eight heterogeneous EEG corpora (712,832 windows; 697,906 with valid labels) to a common 14-channel EPOC+ montage with 63-dimensional spectral features, and we recover the real 1-9 arousal self-assessments for MAHNOB-HCI from session.xml metadata. As a benchmark, Random Forest classifiers are trained on seven source domains and evaluated on the held-out target under both zero-shot and 20%-participant few-shot calibration. The benchmark exposes two concrete methodological pitfalls rather than a performance result: (i) per-class recall shows every target collapsing to a single majority class, and (ii) a within-dataset upper-bound experiment (Table 3) shows that six of eight proxy label sets sit at or below three-class chance even when trained and tested on the same dataset, so the cross-dataset failure is a label-validity problem rather than a transfer-method problem. Across the eight targets (20 seeds, 8,000 evaluation windows per target), zero-shot accuracy averages 36.85% (95% CI 34.40-39.30) and calibrated 43.76% (41.77-45.75), but balanced accuracy stays at 33.01-35.62% (Cohen's kappa <= 0.068), i.e. at chance. The +6.91pp mean change is driven almost entirely by a single target, ds006437 (6.31% -> 60.60%): the median paired change across all 160 seed-pairs is 0.00pp, and after Holm-Bonferroni correction only ds006437 and ds004572 remain significant, the latter with a practically null effect (+0.39pp). Balanced accuracy stays between 33.01% and 35.62% and Cohen's kappa at 0.009 +/- 0.032, i.e. at or barely above three-class chance, while per-class recall shows six of eight targets collapsing to Deep (96.7-100% recall) and two to Light (68.5-99.1%). The collapse persists under SMOTE oversampling, under an EEGNet-v4 deep-learning baseline, and under CORAL and AdaBN feature alignment, which locates the bottleneck in proxy-label validity and class overlap in the feature space rather than in classifier capacity. We position this work as a preliminary methodological study: its contribution is a reproducible eight-dataset alignment pipeline, recovered MAHNOB-HCI arousal self-assessments, a quantitative estimate of split-leakage inflation, and a transparently reported negative result rather than a performance claim.

Inferring "whether a change in the expression of a given gene causally affects the disease state" from observational single-cell transcriptomic data is one of the central problems in single-cell biology. The difficulty lies in confounding: cell state, batch, cell cycle, and the co-expression of other genes may all simultaneously influence the target gene (treatment variable T) and the disease label (outcome variable Y), so that naive correlation analysis cannot distinguish causation from covariation. Double machine learning (DML), via orthogonal scores and cross-fitting, allows machine learning to estimate high-dimensional nuisance functions, thereby addressing the causal inference problem in high-dimensional data. Constrained by computational resources, this study takes a small-sample dataset with p{approx}n (2,120 cells, 1,999 background genes) as the experimental testbed and systematically compares three nuisance function construction strategies under this critical condition; strategies for the n>>p regime are then addressed by theoretical argument. Using systemic lupus erythematosus (SLE) peripheral blood memory B cells (GSE189050, 2,120 cells), we construct a three nuisance function construction strategies x (in-sample / cross-fitting) 2x3 factorial experiment and compare them in terms of resolution, biological plausibility, stability, and deconfounding ability in causal effect estimation. The three strategies are: (S1) direct linear nuisance regression on the high-dimensional background genes without dimensionality reduction; (S2) learning a shared low-dimensional representation with an unsupervised autoencoder, with the treatment and outcome residuals sharing that representation; (S3) fitting the treatment and outcome with two independent deep networks. The results yield a clear three-part pattern with mechanistic meaning: (1) both modes of S1 fail; (2) S2 achieves the optimum in-sample and requires no cross-fitting; (3) S3 is "rescued" under cross-fitting and becomes effective. We explain this pattern starting from the convergence rate condition of DML (Sec.5 Theoretical Foundations) and give the applicability boundaries of the two strategies: S2 (shared unsupervised deep learning) achieves estimation quality comparable to standard cross-fitted DML at a tiny fraction of the computational cost, making it a feasible scheme for large-scale screening; S3 (dual networks + cross-fitting), as the standard DML recipe, can serve as a broad-spectrum reference control for S2 results.

Recurrent neural networks (RNNs) provide a central tool in neuroscience to optimize cognitive tasks and reconstruct dynamical systems. Around these two paradigms, many RNN model variants have been tailored to capture specific computations of the brain. In practice, however, this methodological diversity remains difficult to exploit, because different paradigms and model variants are mostly implemented in idiosyncratic, study-specific code. Here we introduce NeuralRNN, an open-source Python framework that unifies these paradigms and model variants in a general-purpose pipeline. NeuralRNN brings task optimization, dynamical system reconstruction, and biological constraints into a single objective interface. It also provides an integrated and extensible model library spanning generic, structure-constrained, gain-modulated, and gated RNNs. Furthermore, data, models, and training are constructed through a shared automated pipeline and are equipped with a suite of tools for dynamical systems analysis. Together, NeuralRNN lowers the barrier to implementing, comparing, and extending RNN methods in neuroscience.

Subject-independent affect regression from physiological signals remains challenging because emotional responses vary substantially across individuals, widely used datasets provide only coarse trial-level annotations, and heterogeneous physiological modalities may not contribute reliably when treated as if they were interchangeable predictors. We have developed AffectRoute, a protocol-conditioned subject-independent affect regression that is conditioned on the protocol and assigns separate predictive functions to information from the population, electroencephalography (EEG), and peripheral physiological signals (PPS). First, a source-population prior establishes a trial-level affective anchor using only the data from source participants. TrajBridge then combines an EEG representation that is supervised by REFED for participant-specific adjustments with temporal structure obtained from the continuous REFED annotations in order to create a weakly supervised segment-resolved pseudo-trajectory and to establish a frozen trial-level baseline. PhysioRoute next reduces the remaining error by breaking down the residual correction into a source-derived direction, which is estimated from the out-of-fold residuals within the source group, and a channel-specific magnitude derived from the PPS. When evaluated on DEAP and DREAMER using a leave-one-subject-out approach at the participant level, AffectRoute showed consistent step-by-step improvements in both the mean absolute error and the concordance correlation coefficient. A method that relied solely on the source data was clearly worse than PhysioRoute, showing that the final improvement cannot be accounted for by transferable source residual regularity alone. Conventional alternatives to fusing the PPS were also found to be consistently less effective, although analyses at the channel level and with a leave-one-channel-out design showed that the peripheral contributions are axis-dependent yet distributed across channels. These results indicate that structured residual inference is an effective alternative to unrestricted multimodal fusion for subject-independent affect regression.